float a;
float dy = 16;
int k = 0;
int detectorAngles = 9;
int ellipsex;
int ellipsey;
int ellipseflag = 0;
int numvertices = 3;
int inc = 0;
int numRays;
float[] polydata = new float[2*numvertices];
int[][] testRay;
Reconstruction[] backProjection = new Reconstruction[detectorAngles];

void setup() 
{
  size(1536, 512);
  stroke(255);
  smooth();
  numRays = height;
  a = height;
  testRay = new int[detectorAngles][numRays+1];
  background(0);
}
void mousePressed() {
       ellipsex = mouseX;
       ellipsey = mouseY; 
       ellipseflag = 1;
     }

void Polygon(){
    beginShape();
    for(int m = 0; m < numvertices; m=m+1) {
      vertex(ellipsex/2*cos(2*PI*m/numvertices),ellipsey/2*sin(2*m*PI/numvertices));
      polydata[2*m] = ellipsex/2*cos(2*PI*m/numvertices);
      polydata[2*m+1] = ellipsey/2*sin(2*PI*m/numvertices);
    }
    endShape(CLOSE);
}

void draw() 
{
  float angle;
  float xtest, ytest;
  int insidetest = 0;
  int jj;
  int numRays = 2*int(height) / int(dy);

  fill(0);
  rect(0,0,width/3,height);
  line(2*width/3,0,2*width/3,height);

  if (ellipseflag==0){
    print("Use the mouse to select the ellipse size, referenced to top left corner\n");
  }
  if (mousePressed == true) {
      Polygon();
  }
  if (ellipseflag==1){   
      if (k == detectorAngles)
      {
        k = 0;
        noLoop();
      }
    if(k<detectorAngles){
      translate(width/6, height/2);
      fill(255);
      Polygon();
      angle = k*PI/detectorAngles;
      rotate(angle);
      line(-width/12, a/4, width/12, a/4);  
      for (int g = 0; g < 64; g=g+1){
        xtest = -width/12 + g/64.0 * width/6;
        ytest = a/4;      
        insidetest=insidetest + insidePolygon(xtest,ytest,polydata,angle);        
      }
    }
    testRay[k][inc] = insidetest;
    a = a - dy;
    inc = inc + 1;
    if (a < -height) { 
      a = height; 
      inc = 0;
      if(k<=detectorAngles){ 
        backProjection[k] = new Reconstruction();
        rotate(-k*PI/detectorAngles);
        translate(width/3, 0*height/2);
        backProjection[k].display(0);
        rotate(-k*PI/detectorAngles);
        translate(width/3, 0*height/2);
        backProjection[k].display(1);
      }
      k = (k + 1);
    }
  }
}

int insidePolygon(float x,float y, float poly[],float rotangle) {
  int pnum=poly.length/2;
  int id1,id2;
  int i, j, c = 0;
  float[] polyr = new float[2*pnum];
  
  for(int n = 0; n < pnum; n=n+1){
    polyr[2*n] = poly[2*n]*cos(-rotangle) - poly[2*n+1]*sin(-rotangle);
    polyr[2*n+1] = poly[2*n]*sin(-rotangle) + poly[2*n+1]*cos(-rotangle);
  }    

  for (i = 0, j = pnum-1; i < pnum; j = i++) {
    id1=i*2;
    id2=j*2;
    if ((((polyr[id1+1] <= y) && (y < polyr[id2+1])) || ((polyr[id2+1] <= y) && (y < polyr[id1+1]))) &&
	     (x < (polyr[id2] - polyr[id1]) * (y - polyr[id1+1]) / (polyr[id2+1] - polyr[id1+1]) + polyr[id1]))
	     c = (c+1)%2;
  }
  if(c==1){ 
    stroke(255,0,0);
    fill(0,0,0);
    ellipse(x,y,4,4);
    stroke(255);
    return 2;
  } else {
    return 0;
  }
}

class Reconstruction
{
  int numRays = 2*int(height) / int(dy);
  int reconType;
  int cl;

  Reconstruction(){}
  
  void display(int type){
    float Opacity = 0.4*255;    
    reconType = type;
    rotate(k*PI/detectorAngles);
    for (int g = 0; g < numRays; g=g+1){
      if (reconType==1){
        if(testRay[k][g]>1.0){
          stroke(2*testRay[k][g],Opacity);
        } else {
          stroke(25,Opacity);
        }
      } else {
        if(testRay[k][g]>1.0){
          stroke(255,Opacity);
        } else {
          stroke(50,Opacity);
        }
      }
      line(-width/12, (a-g*dy)/4, width/12, (a-g*dy)/4);            
    }
  stroke(255);
  }
}