import java.io.*;
import java.util.*;

import javax.imageio.ImageIO;
import javax.imageio.ImageReadParam;
import javax.imageio.ImageReader;

import edu.rit.pj.*;
import edu.rit.matrix.ByteArray;
import edu.rit.matrix.ByteMatrix;
import edu.rit.matrix.GrayImage;

public class EdgeLineSmp {
    private static GrayImage readImage(String infile) throws IOException {
        String suffix = infile.substring(infile.lastIndexOf('.') + 1);
        Iterator i = ImageIO.getImageReadersBySuffix(suffix);
        if(!i.hasNext()) throw new IOException("no reader for that image type");
        ImageReader reader = (ImageReader) i.next();        
        reader.setInput(ImageIO.createImageInputStream(new File(infile)),true);        
        int w = reader.getWidth (0);
        int h = reader.getHeight (0);
        ByteMatrix matrix = new ByteMatrix (w, h);
        matrix.allocate();        
        GrayImage image = new GrayImage(matrix);        
        ImageReadParam param = reader.getDefaultReadParam();
        param.setDestination(image);
        reader.read(0, param);
        return image;
    }
    
    public static void main(String[] args) throws Exception {
        Comm.init (args);
        
        if(args.length != 4) throw new IllegalArgumentException("wrong number of arguments");
        final long t = -System.currentTimeMillis();
        GrayImage image = readImage(args[0]);
        final double dtheta, dr, thr;
        try {
            dtheta = Double.parseDouble(args[1]);
            dr = Double.parseDouble(args[2]);
            thr = Double.parseDouble(args[3]);
        } catch(NumberFormatException e) {
            throw new IllegalArgumentException("invalid double value");
        }
        
        int w = image.getWidth();
        int h = image.getHeight();
        double[] pixels = new double[w*h];
        for(int i=0;i<h;i++) for(int j=0;j<w;j++) pixels[i*w+j] = image.getPixel(i,j);
        
        boolean[] isEdgePixel = new boolean[(w-1)*(h-1)];
        for(int i=0;i<h-1;i++)
            for(int j=0;j<w-1;j++)
                isEdgePixel[i*(w-1)+j] = Math.abs(pixels[i*w+j] - pixels[(i+1)*w+(j+1)]) + Math.abs(pixels[(i+1)*w+j] - pixels[i*w+(j+1)]) != 0.0;
                
        int edgePixelCount_ = 0;
        for(int i=0;i<isEdgePixel.length;i++) if(isEdgePixel[i]) edgePixelCount_++;
        
        final int[] edgePixelsX = new int[edgePixelCount_];
        final int[] edgePixelsY = new int[edgePixelCount_];
        edgePixelCount_ = 0;
        
        for(int i=0;i<h-1;i++) {
            for(int j=0;j<w-1;j++) {
                if(isEdgePixel[i*(w-1)+j]) {
                    edgePixelsX[edgePixelCount_] = j;
                    edgePixelsY[edgePixelCount_] = i;
                    edgePixelCount_++;
                }
            }
        }
        final int edgePixelCount = edgePixelCount_;
        if(edgePixelCount != edgePixelsX.length) throw new Error("should never happen");
        
        final int numThetaBins = (int) Math.ceil(Math.PI * 2 / dtheta);
        final int numRBins = (int) Math.ceil(Math.sqrt((w-1)*(w-1) + (h-1)*(h-1)) / dr);
        final int numBins = numThetaBins * numRBins;
        
        final double[] thetaSum = new double[numBins];
        final double[] rSum = new double[numBins];
        final int[] count = new int[numBins];
        
        new ParallelTeam().execute (new ParallelRegion() {
            public void run() throws Exception {
                final double[] lthetaSum = new double[numBins];
                final double[] lrSum = new double[numBins];
                final int[] lcount = new int[numBins];
                
                execute(0,edgePixelCount-1,new IntegerForLoop() {
                    public Schedule<Integer> schedule() { return IntegerSchedule.guided(256); }
                    public void run(int first, int last) {
                        for(int i=first;i<=last;i++) {
                            int x1 = edgePixelsX[i];
                            int y1 = edgePixelsY[i];
                            for(int j=i+1;j<edgePixelCount;j++) {
                                int x2 = edgePixelsX[j];
                                int y2 = edgePixelsY[j];
                                int dx = x2 - x1;
                                int dy = y2 - y1;
                                double theta,r;
                                if(dx == 0) {
                                    theta = 0.0;
                                    r = x1 + 0.5;
                                } else if(dy ==  0) {
                                    theta = Math.PI / 2;
                                    r = y1 + 0.5;
                                } else {
                                    double m = (double)dy / (double) dx;
                                    double b = (y1 + 0.5) - m*(x1 + 0.5);
                                    theta = Math.atan(-1/m);
                                    r = Math.sin(theta) * b;
                                    if(r < 0) { r = -r; theta += Math.PI; }
                                    else if(theta < 0) theta += Math.PI*2;
                                }
                                
                                int bin = ((int)(theta / dtheta)) * numRBins + ((int)(r/dr));
                                lthetaSum[bin] += theta;
                                lrSum[bin] += r;
                                lcount[bin] += 1;
                            }
                        }
                    }
                });
                critical(new ParallelSection() {
                    public void run() {
                        for(int i=0;i<numBins;i++) {
                            thetaSum[i] += lthetaSum[i];
                            rSum[i] += lrSum[i];
                            count[i] += lcount[i];
                        }
                    }
                });
            }
        });
        
        
        int largest = 0;
        for(int i=0;i<numBins;i++) {
            int c = count[i];
            if(c > largest) largest = c;
        }
        int large = (int) Math.floor(thr * largest);
        if(large != 0) for(int i=0;i<numBins;i++) {
            int c = count[i];
            if(c >= large) {
                double theta = thetaSum[i] / c;
                double r = rSum[i] / c;
                if(theta > Math.PI) theta = theta - Math.PI*2;
                System.out.println(theta + "\t" + r);
            }
        }
        

        System.err.println(t+System.currentTimeMillis());
    }
}