#include <iostream>
#include <fstream>
#include <math.h>
#include <ctype.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>

static const double SMALL = 0.0000001;
static const int SMALLEST_SET = 25;

static const double hueWeight = 0.25;
static const double satWeight = 0.15;
static const double valWeight = 1.0 - hueWeight - satWeight;

    /*
    ** Something must be minRatio as tall as it is wide to
    ** be considered.
    */
static const double smallCutoff = 60;
static const double smallRatio = 0.65;
static const double bigRatio = 0.42;



static bool
ReadPPMHeader( istream& in, unsigned int* w, unsigned int* h )
{
    bool gotP = false;
    bool gotSix = false;
    bool inWidth = false;
    bool gotWidth = false;
    bool inHeight = false;
    bool gotHeight = false;
    bool inMax = false;
    bool gotMax = false;

    *w = 0;
    *h = 0;

    in.unsetf( ios::skipws );

    while (!gotMax) {
	char c;
	in >> c;

	if ( !in.good() ) {
	    in.setf( ios::skipws );
	    return false;
	}

	if ( c == '#' ) {
	    while ( c != '\n' ) {
		in >> c;
		if ( !in.good() ) {
		    in.setf( ios::skipws );
		    return false;
		}
	    }
	} else if ( !gotP ) {
	    if ( c == 'P' ) {
		gotP = true;
	    }
	} else if ( !gotSix ) {
	    if ( c == '6' ) {
		gotSix = true;
	    }
	} else if ( isdigit( c ) ) {
	    if ( !gotWidth ) {
		inWidth = true;
		*w = *w * 10 + c - '0';
	    } else if ( !gotHeight ) {
		inHeight = true;
		*h = *h * 10 + c - '0';
	    } else if ( !gotMax ) {
		inMax = true;
	    }
	} else {
	    if ( inWidth ) {
		gotWidth = true;
		inWidth = false;
	    } else if ( inHeight ) {
		gotHeight = true;
		inHeight = false;
	    } else if ( inMax ) {
		gotMax = true;
		inMax = false;
	    }
	}
    }

    return true;
}

static inline unsigned char
clip( double v )
{
    return (unsigned char)( 
	    ((v) < 0.0)
		? 0.0 : ((v) > 255.0)
			    ? 255.0 : v
	);
}


static void
ConvertToHSV(
	const unsigned char* im,
	unsigned int width, unsigned int height,
	double* hh, double* ss, double* vv
    )
{
#ifndef NDEBUG
    const unsigned char* saveIm = im;
    const double* saveHH = hh;
    const double* saveSS = ss;
    const double* saveVV = vv;
#endif

    for ( unsigned int jj=0; jj < height; ++jj ){
	for ( unsigned int ii=0; ii < width; ++ii ) {
	    double rr = ( (double)*im++ ) / 255.0;
	    double gg = ( (double)*im++ ) / 255.0;
	    double bb = ( (double)*im++ ) / 255.0;

	    double mx = rr;
	    double mn = rr;
	    unsigned int mxSep = 0;

	    if ( gg > mx ) {
		mx = gg;
		mxSep = 1;
	    } else if ( gg < mn ) {
		mn = gg;
	    }

	    if ( bb > mx ) {
		mx = bb;
		mxSep = 2;
	    } else if ( bb < mn ) {
		mn = bb;
	    }


	    double delta = mx - mn;

	    *vv = mx;

	    if ( mx > SMALL ) {
		*ss = delta / mx;
	    } else {
		*ss = 0.0;
	    }

	    if ( *ss > SMALL ) {
		switch ( mxSep ) {
		case 0:
		    *hh = ( gg - bb ) / delta;
		    break;
		case 1:
		    *hh = 2.0 + ( bb - rr ) / delta;
		    break;
		case 2:
		    *hh = 4.0 + ( rr - gg ) / delta;
		    break;
		}

		*hh *= M_PI / 3.0 ;

		if ( *hh < 0.0 ) {
		    *hh += 2.0 * M_PI;
		}

	    } else {
		*hh = 0.0;
	    }

	    ++hh;
	    ++ss;
	    ++vv;
	}
    }

    assert( ( im - saveIm ) == 3 * width * height );
    assert( ( hh - saveHH ) == width * height );
    assert( ( ss - saveSS ) == width * height );
    assert( ( vv - saveVV ) == width * height );
}

static void
ConvertFromHSV(
	const double* fhh, const double* fss, const double* fvv,
	unsigned int width, unsigned int height,
	unsigned char* im
    )
{
#ifndef NDEBUG
    const double* saveHH = fhh;
    const double* saveSS = fss;
    const double* saveVV = fvv;
    const unsigned char* saveIm = im;
#endif

    for ( unsigned int jj=0; jj < height; ++jj ){
	for ( unsigned int ii=0; ii < width; ++ii ) {
	    double hh = *fhh++;
	    double ss = *fss++;
	    double vv = *fvv++;

		/*
		** make sure all of the values are in
		** the range we want them to be.
		*/
	    while ( hh > 5.0 * M_PI / 3.0 ) {
		hh -= 2.0 * M_PI;
	    }

	    while ( hh < -M_PI / 3.0 - SMALL ) {
		hh += 2.0 * M_PI;
	    }

	    hh /= M_PI / 3.0;

	    if ( ss > 1.0 - SMALL ) {
		ss = 1.0;
	    } else if ( ss < SMALL ) {
		ss = 0.0;
	    }

	    if ( vv > 1.0 - SMALL ) {
		vv = 1.0;
	    } else if ( vv < SMALL ) {
		vv = 0.0;
	    }

		/*
		** Determine rr, gg, and bb
		*/
	    double rr;
	    double gg;
	    double bb;

	    if ( ss > SMALL ) {
		double delta = ss * vv;

		if ( hh <= 2.0 ) {
		    double diff = hh * delta;

		    rr = vv;

		    if ( diff < SMALL ) {
			gg = vv - delta;
			bb = gg - diff;
		    } else {
			bb = vv - delta;
			gg = diff + bb;
		    }

		} else if ( hh <= 4.0 ) {
		    double diff = ( hh - 2.0 )  * delta;

		    gg = vv;

		    if ( diff < SMALL ) {
			bb = vv - delta;
			rr = bb - diff;
		    } else {
			rr = vv - delta;
			bb = diff + rr;
		    }

		} else {
		    double diff = ( hh - 4.0 )  * delta;

		    bb = vv;

		    if ( diff < SMALL ) {
			rr = vv - delta;
			gg = rr - diff;
		    } else {
			gg = vv - delta;
			rr = diff + gg;
		    }
		}
	    } else {
		rr = gg = bb = vv;
	    }

		/*
		** Output pixels
		*/
	    *im++ = clip( rr * 255.0 );
	    *im++ = clip( gg * 255.0 );
	    *im++ = clip( bb * 255.0 );
	}
    }

    assert( ( im - saveIm ) == 3 * width * height );
    assert( ( fhh - saveHH ) == width * height );
    assert( ( fss - saveSS ) == width * height );
    assert( ( fvv - saveVV ) == width * height );
}

static void
ConvertFromH(
	const double* rr,
	unsigned int width, unsigned int height,
	unsigned char* im
    )
{
#ifndef NDEBUG
    const double* saveRR = rr;
    const unsigned char* saveIm = im;
#endif

    for ( unsigned int jj=0; jj < height; ++jj ) {
	for ( unsigned int ii=0; ii < width; ++ii ) {
	    double score = *rr++;

	    *im++ = clip( score * 255.0 );
	    *im++ = clip( score * 255.0 );
	    *im++ = clip( score * 255.0 );
	}
    }

    assert( ( rr - saveRR ) == width * height );
    assert( ( im - saveIm ) == 3 * width * height );
}


static void
ScoreDifferences(
	const double* b0, const double* b1, const double* b2,
	const double* f0, const double* f1, const double* f2,
	unsigned int width, unsigned int height,
	double* rr
    )
{
    for ( unsigned int jj=0; jj < height; ++jj ) {
	for ( unsigned int ii=0; ii < width; ++ii ) {
	    unsigned int index = ii + jj * width;

	    assert( index < width * height );

	    double hueDiff = fabs( f0[ index ] - b0[ index ] );
	    double satDiff = fabs( f1[ index ] - b1[ index ] );
	    double valDiff = fabs( f2[ index ] - b2[ index ] );

	    while ( hueDiff >= 2.0 * M_PI ){
		hueDiff -= 2.0 * M_PI;
	    }

	    if ( hueDiff > M_PI ){
		hueDiff = 2.0 * M_PI - hueDiff;
	    }

	    hueDiff /= M_PI;

	    double score = hueWeight * hueDiff
			+ satWeight * satDiff
			+ valWeight * valDiff;

	    rr[ index ] = pow( score, 0.55 );
	}
    }
}

typedef struct PixelTree {
    unsigned int index;
    double val;
    unsigned int cntr;
    struct PixelTree* less;
    struct PixelTree* more;
} PixelTree;

static inline unsigned int
AddPixel( PixelTree* array, unsigned int cntr, double val )
{
    array[ cntr ].index = cntr;
    array[ cntr ].val = val;
    array[ cntr ].cntr = 1;
    array[ cntr ].less = 0;
    array[ cntr ].more = 0;

    if ( cntr > 0 ) {
	PixelTree** ptrPtr = 0;
	PixelTree* ptr = &array[ 0 ];

	while ( ptr != 0 ) {
	    assert( ptr->index < cntr );

	    ++ptr->cntr;

	    if ( val < ptr->val ) {
		ptrPtr = &ptr->less;
	    } else {
		ptrPtr = &ptr->more;
	    }

	    ptr = *ptrPtr;
	}

	*ptrPtr = &array[ cntr ];
    }

    return cntr+1;
}

static inline double
GetMedian( PixelTree* ptr )
{
    unsigned int targetCnt = ptr->cntr / 2 + 1;

    while ( ptr->cntr > 1 ) {
	if ( ptr->less != 0 ) {
	    if ( ptr->less->cntr >= targetCnt ) {
		ptr = ptr->less;
	    } else if ( ptr->less->cntr+1 == targetCnt ) {
		break;
	    } else {
		targetCnt -= ptr->less->cntr+1;
		ptr = ptr->more;
	    }
	} else {
	    if ( targetCnt == 1 ) {
		break;
	    } else {
		--targetCnt;
		ptr = ptr->more;
	    }
	}
    }

    return ptr->val;
}

static void
MedianFilter(
	const double* in, unsigned int width, unsigned int height,
	unsigned int rad, double* out
    )
{
    unsigned int rad2 = rad * rad;

    PixelTree* tmp = new PixelTree[ rad2 * rad2 ];

    for ( unsigned int jj=0; jj < height; ++jj ) {
	for ( unsigned int ii=0; ii < width; ++ii ) {
	    unsigned int cntr = 0;

	    for ( int yy = -(int)rad; yy <= (int)rad; ++yy ) {
		int ycoord = yy + (int)jj;
		int xrad2 = rad2 - yy * yy;

		if ( ycoord < 0 || ycoord >= height ){
		    continue;
		}

		for ( int xx = -(int)rad; xx <= (int)rad; ++xx ) {
		    int xcoord = xx + (int)ii;

		    if ( xcoord < 0 || xcoord >= width ){
			continue;
		    }

		    if ( xx * xx > xrad2 ) {
			continue;
		    }

		    assert( xcoord + ycoord * width < width * height );
		    cntr = AddPixel( tmp, cntr, in[ xcoord + ycoord * width ] );
		}
	    }

	    double res = in[ ii + jj * width ];

	    if ( cntr > 0 ) {
		res = GetMedian( tmp );
	    }

	    *out++ = res;
	}
    }

    delete[] tmp;
}

static void
RescaleValues(
	const double* in,
	unsigned int width, unsigned int height,
	double* out
    )
{
    for ( unsigned int jj=0; jj < height; ++jj ) {
	for ( unsigned int ii=0; ii < width; ++ii ) {
	    double val = *in++ * 255.0;

	    if ( val < 115.5 ) {
		val = 0.0;
	    } else if ( val < 127.5 ) {
		val = ( val - 115.0 ) * ( 173.0-80.0 ) / ( 127.0 - 115.0 )
			+ 80.0;
	    } else if ( val < 142.5 ) {
		val = ( val - 127.0 ) * ( 228.0-173.0 ) / ( 142.0-127.0 )
			+ 173.0;
	    } else {
		val = ( val - 142.0 ) * ( 255.0-228.0 ) / ( 255.0-142.0 )
			+ 228.0;
	    }

	    *out++ = val / 255.0;
	}
    }
}

static void
ApplyMask(
	const double* mask,
	unsigned int width, unsigned int height,
	unsigned char* im
    )
{
    for ( unsigned int jj=0; jj < height; ++jj ) {
	for ( unsigned int ii=0; ii < width; ++ii ) {
	    if ( *mask++ > SMALL ) {
		im += 3;
	    } else {
		*im++ |= 0x80;
		*im++ |= 0xF0;
		*im++ |= 0x80;
	    }
	}
    }
}

static inline unsigned int
SetOf( int* sets, unsigned int index )
{
    unsigned int parent = index;

    while ( sets[ parent ] >= 0 ) {
	parent = sets[ parent ];
    }

	/*
	** Path compression
	*/
    while ( sets[ index ] >= 0 ) {
	unsigned int saveIndex = index;
	index = sets[ index ];
	sets[ saveIndex ] = parent;
    }

    return parent;
}

static inline unsigned int
SetCount( const int* sets, unsigned int index )
{
    unsigned int parent = SetOf( (int*)sets, index );
    return (unsigned int)-sets[ parent ];
}

static inline void
JoinSets( int* sets, unsigned int a, unsigned int b )
{
    a = SetOf( sets, a );
    b = SetOf( sets, b );

    if ( a < b ) {
	sets[ a ] += sets[ b ];
	sets[ b ] = a;
    } else if ( b < a ) {
	sets[ b ] += sets[ a ];
	sets[ a ] = b;
    }
}

static unsigned int
SegmentImage(
	const double* r0,
	unsigned int width,
	unsigned int height,
	int* sets
    )
{
    for ( unsigned int ii=0; ii < width * height; ++ii ) {
	sets[ii] = -1;
    }

    const double* rCur = r0;
    unsigned int index = 0;

    for ( unsigned int jj=0; jj < height; ++jj ) {
	const double* rPrev = 0;

	for ( unsigned int ii=0; ii < width; ++ii ) {
	    if ( *rCur > SMALL ) {
		unsigned int ss;

		if ( ii > 0 ) {
		    if ( ( ss=SetOf( sets, index-1 ) ) != 0 ) {
			JoinSets( sets, ss, index );
		    }

		    if ( jj > 0 ) {
			if ( ( ss=SetOf( sets, index-1-width ) ) != 0 ) {
			    JoinSets( sets, ss, index );
			}
		    }
		}

		if ( jj > 0 ) {
		    if ( ( ss=SetOf( sets, index+0-width ) ) != 0 ) {
			JoinSets( sets, ss, index );
		    }
		}

		if ( ii+1 < width && jj > 0 ) {
		    if ( ( ss=SetOf( sets, index+1-width ) ) != 0 ) {
			JoinSets( sets, ss, index );
		    }
		}
	    } else {
		JoinSets( sets, 0, index );
	    }

	    ++rCur;
	    ++index;
	}
    }

    unsigned int uniqueSets = 0;

    for ( unsigned int ii=0; ii < width * height; ++ii ) {
	if ( sets[ii] < 0 ) {
		/*
		** Filter out really small bits
		** count bigger bits
		*/
	    if ( ii == 0 || sets[ii] < -SMALLEST_SET ) {
		sets[ii] = - (int) ++uniqueSets;
	    } else {
		sets[ii] = 0;
	    }
	}
    }

    return uniqueSets;
}

static void
PlotSets(
	const int* sets, unsigned int uniqueSets, unsigned int totalSets,
	double* r1
    )
{
    double scale = 0.0;

    if ( uniqueSets > 1 ) {
	scale = 1.0 / (double)( uniqueSets - 1 );
    }

    for ( unsigned int ii=0; ii < totalSets; ++ii ) {
	*r1++ = ( (double) SetCount( sets, ii ) - 1.0 ) * scale;
    }
}

typedef struct {
    unsigned int cntr;
    unsigned int lastY;
    unsigned int x;
    unsigned int y;
    unsigned int minX;
    unsigned int maxX;
    unsigned int minY;
    unsigned int maxY;
} CentroidCalc;

static bool
RatioTest( CentroidCalc* cc )
{
    double dy = (double)( cc->maxY - cc->minY );
    double dx = (double)( cc->maxX - cc->minX );

    if ( cc->cntr < smallCutoff ) {
	return ( dy > smallRatio * dx );
    } else {
	return ( dy > bigRatio * dx );
    }
}

static void
CentroidSpecificSetWithDivider(
	int* sets, unsigned int set,
	unsigned int width, unsigned int height,
	unsigned int divx, unsigned int cutoff
    )
{
    CentroidCalc calc[2];

    memset( calc, 0, 2 * sizeof( CentroidCalc ) );

    calc[0].minX = width;
    calc[0].minY = height;
    calc[1].minX = width;
    calc[1].minY = height;

    for ( unsigned int jj=0; jj < height; ++jj ) {
	for ( unsigned int ii=0; ii < width; ++ii ) {
	    unsigned int ss = SetCount( sets, ii + jj * width ) - 1;

	    if ( ss == set ) {
		unsigned int kk = 0;

		if ( ii > divx ) {
		    kk = 1;
		}

		if ( calc[kk].lastY == 0 ) {
		    calc[kk].lastY = jj + cutoff;
		}

		if ( ii < calc[kk].minX ) {
		    calc[kk].minX = ii;
		} else if ( ii > calc[ss].maxX ) {
		    calc[kk].maxX = ii;
		}

		if ( jj < calc[kk].minY ) {
		    calc[kk].minY = jj;
		} else if ( jj > calc[ss].maxY ) {
		    calc[kk].maxY = jj;
		}

		if ( jj < calc[kk].lastY ) {
		    ++calc[kk].cntr;
		    calc[kk].x += ii;
		    calc[kk].y += jj;
		}
	    }
	}
    }

    for ( unsigned int ii=0; ii < 2; ++ii ) {
	if ( calc[ii].cntr > 0 ) {
	    unsigned int xx = calc[ii].x / calc[ii].cntr;
	    unsigned int yy = calc[ii].y / calc[ii].cntr;

	    unsigned int ss = SetCount( sets, xx + yy * width ) - 1;

	    if ( ss == set ) {
		if ( yy > height/16 && RatioTest( &calc[ii] ) ) {
		    cout << xx << " " << yy << endl;
		    sets[ xx + yy * width ] = 0;
		} else {
		    // JoinSets( sets, 0, xx + yy * width );
		}
	    }
	}
    }
}

static void
FindCentroidsOfSets(
	int* sets, unsigned int uniqueSets,
	unsigned int width, unsigned int height
    )
{
    CentroidCalc* calc = new CentroidCalc[ uniqueSets ];

    memset( calc, 0, uniqueSets * sizeof( CentroidCalc ) );

    for ( unsigned int ii=0; ii < uniqueSets; ++ii ) {
	calc[ii].minX = width;
	calc[ii].minY = height;
    }

    unsigned int cutoff = height/6;

    for ( unsigned int jj=0; jj < height; ++jj ) {
	for ( unsigned int ii=0; ii < width; ++ii ) {
	    unsigned int ss = SetCount( sets, ii + jj * width ) - 1;

	    if ( calc[ss].lastY == 0 ) {
		calc[ss].lastY = jj + cutoff;
	    }

	    if ( ii < calc[ss].minX ) {
		calc[ss].minX = ii;
	    } else if ( ii > calc[ss].maxX ) {
		calc[ss].maxX = ii;
	    }

	    if ( jj < calc[ss].minY ) {
		calc[ss].minY = jj;
	    } else if ( jj > calc[ss].maxY ) {
		calc[ss].maxY = jj;
	    }

	    if ( jj < calc[ss].lastY ) {
		++calc[ss].cntr;
		calc[ss].x += ii;
		calc[ss].y += jj;
	    }
	}
    }

	/*
	** Skip the zero-th set.
	*/
    for ( unsigned int ii=1; ii < uniqueSets; ++ii ) {
	if ( calc[ii].cntr > 0 ) {
	    unsigned int xx = calc[ii].x / calc[ii].cntr;
	    unsigned int yy = calc[ii].y / calc[ii].cntr;

	    unsigned int ss = SetCount( sets, xx + yy * width ) - 1;

	    if ( ss == ii ) {
		if ( yy > height/16 && RatioTest( &calc[ii] ) ) {
		    cout << xx << " " << yy << endl;
		    sets[ xx + yy * width ] = 0;
		} else {
		    // JoinSets( sets, 0, xx + yy * width );
		}
	    } else {
		CentroidSpecificSetWithDivider(
			sets, ii, width, height, xx, cutoff
		    );
	    }
	}
    }

    delete[] calc;
}

int
main( int argc, const char* argv[] )
{
    unsigned int width;
    unsigned int height;

    ifstream baseFile( argv[1], ios::in );

    if ( !ReadPPMHeader( baseFile, &width, &height ) ) {
	cerr << "Cannot read PPM header" << endl;
	return __LINE__;
    }

    unsigned int fwidth;
    unsigned int fheight;

    ifstream findFile( argv[2], ios::in );

    if ( !ReadPPMHeader( findFile, &fwidth, &fheight ) ) {
	cerr << "Cannot read PPM header" << endl;
	return __LINE__;
    }

    if ( fwidth != width || fheight != height ) {
	cerr << "Images must be the same dimensions" << endl;
	return __LINE__;
    }

    unsigned char* base = new unsigned char[ 3*width*height ];
    baseFile.read( base, 3*width*height );

    unsigned char* find = new unsigned char[ 3*width*height ];
    findFile.read( find, 3*width*height );

    double* b0 = new double[ width*height ];
    double* b1 = new double[ width*height ];
    double* b2 = new double[ width*height ];

    double* f0 = new double[ width*height ];
    double* f1 = new double[ width*height ];
    double* f2 = new double[ width*height ];

    double* r0 = new double[ width*height ];
    double* r1 = new double[ width*height ];

    int* sets = new int[ width * height ];

    ConvertToHSV( base, width, height, b0, b1, b2 );
    ConvertToHSV( find, width, height, f0, f1, f2 );

    ScoreDifferences( b0, b1, b2, f0, f1, f2, width, height, r0 );

    MedianFilter( r0, width, height, 6, r1 );

    RescaleValues( r1, width, height, r0 );

    // ApplyMask( r0, width, height, find );

    unsigned int uniqueSets = SegmentImage( r0, width, height, sets );

    FindCentroidsOfSets( sets, uniqueSets, width, height );

#ifndef NDEBUG
    PlotSets( sets, uniqueSets, width * height, r0 );
    ApplyMask( r0, width, height, find );
    // ConvertFromH( r0, width, height, find );

    ofstream dbg( "debug.ppm", ios::out | ios::binary );
    dbg << "P6" << endl;
    dbg << "# CREATOR: people find (debug)" << endl;
    dbg << width << " " << height << endl;
    dbg << 255 << endl;
    dbg.write( find, 3*width*height );
    dbg.close();
#endif

    delete[] sets;

    delete[] r1;
    delete[] r0;

    delete[] f2;
    delete[] f1;
    delete[] f0;

    delete[] b2;
    delete[] b1;
    delete[] b0;

    delete[] find;

    return 0;
}