#include <bits/stdc++.h>
// <graph/ufds.hpp>
#include <optional>
#include <utility>
#include <vector>
template <class Container = std::vector<int>, bool path_compression = true>
struct UFDS {
  Container a;
  UFDS(int n) : a(n, -1) {}
  int find(int i) {
    if constexpr (path_compression) {
      return a[i] < 0 ? i : a[i] = find(a[i]);
    } else {
      while (a[i] >= 0) i = a[i];
      return i;
    }
  }
  std::optional<std::pair<int, int>> unite(int i, int j) {
    i = find(i), j = find(j);
    if (i == j) return std::nullopt;
    if (-a[i] < -a[j]) std::swap(i, j);
    a[i] = a[i] + a[j];
    a[j] = i;
    return std::make_pair(i, j);
  }
};
// </graph/ufds.hpp>

#define rep(i, a, b) for (int i = (a); i < (b); ++i)

using std::ostream;
using std::vector;
using std::tie;

template<class T>
auto sz(T &&x) -> int {
  return static_cast<int>(std::size(x));
}

template <class T> int sgn(T x) { return (x > 0) - (x < 0); }
template<class T>
struct Point {
	typedef Point P;
	T x, y;
	explicit Point(T x_=0, T y_=0) : x(x_), y(y_) {}
	bool operator<(P p) const { return tie(x,y) < tie(p.x,p.y); }
	bool operator==(P p) const { return tie(x,y)==tie(p.x,p.y); }
	P operator+(P p) const { return P(x+p.x, y+p.y); }
	P operator-(P p) const { return P(x-p.x, y-p.y); }
	P operator*(T d) const { return P(x*d, y*d); }
	P operator/(T d) const { return P(x/d, y/d); }
	T dot(P p) const { return x*p.x + y*p.y; }
	T cross(P p) const { return x*p.y - y*p.x; }
	T cross(P a, P b) const { return (a-*this).cross(b-*this); }
	T dist2() const { return x*x + y*y; }
	double dist() const { return sqrt((double)dist2()); }
	// angle to x-axis in interval [-pi, pi]
	double angle() const { return atan2(y, x); }
	P unit() const { return *this/dist(); } // makes dist()=1
	P perp() const { return P(-y, x); } // rotates +90 degrees
	P normal() const { return perp().unit(); }
	// returns point rotated 'a' radians ccw around the origin
	P rotate(double a) const {
		return P(x*cos(a)-y*sin(a),x*sin(a)+y*cos(a)); }
	friend ostream& operator<<(ostream& os, P p) {
		return os << "(" << p.x << "," << p.y << ")"; }
};


template<class T> struct Point3D {
	typedef Point3D P;
	typedef const P& R;
	T x, y, z;
	explicit Point3D(T x_=0, T y_=0, T z_=0) : x(x_), y(y_), z(z_) {}
	bool operator<(R p) const {
		return tie(x, y, z) < tie(p.x, p.y, p.z); }
	bool operator==(R p) const {
		return tie(x, y, z) == tie(p.x, p.y, p.z); }
	P operator+(R p) const { return P(x+p.x, y+p.y, z+p.z); }
	P operator-(R p) const { return P(x-p.x, y-p.y, z-p.z); }
	P operator*(T d) const { return P(x*d, y*d, z*d); }
	P operator/(T d) const { return P(x/d, y/d, z/d); }
	T dot(R p) const { return x*p.x + y*p.y + z*p.z; }
	P cross(R p) const {
		return P(y*p.z - z*p.y, z*p.x - x*p.z, x*p.y - y*p.x);
	}
	T dist2() const { return x*x + y*y + z*z; }
	double dist() const { return sqrt((double)dist2()); }
	//Azimuthal angle (longitude) to x-axis in interval [-pi, pi]
	double phi() const { return atan2(y, x); } 
	//Zenith angle (latitude) to the z-axis in interval [0, pi]
	double theta() const { return atan2(sqrt(x*x+y*y),z); }
	P unit() const { return *this/(T)dist(); } //makes dist()=1
	//returns unit vector normal to *this and p
	P normal(P p) const { return cross(p).unit(); }
	//returns point rotated 'angle' radians ccw around axis
	P rotate(double angle, P axis) const {
		double s = sin(angle), c = cos(angle); P u = axis.unit();
		return u*dot(u)*(1-c) + (*this)*c - cross(u)*s;
	}
};

typedef Point3D<double> P3;

struct PR {
	void ins(int x) { (a == -1 ? a : b) = x; }
	void rem(int x) { (a == x ? a : b) = -1; }
	int cnt() { return (a != -1) + (b != -1); }
	int a, b;
};

struct F { P3 q; int a, b, c; };

vector<F> hull3d(const vector<P3>& A) {
	assert(sz(A) >= 4);
	vector<vector<PR>> E(sz(A), vector<PR>(sz(A), {-1, -1}));
#define E(x,y) E[f.x][f.y]
	vector<F> FS;
	auto mf = [&](int i, int j, int k, int l) {
		P3 q = (A[j] - A[i]).cross((A[k] - A[i]));
		if (q.dot(A[l]) > q.dot(A[i]))
			q = q * -1;
		F f{q, i, j, k};
		E(a,b).ins(k); E(a,c).ins(j); E(b,c).ins(i);
		FS.push_back(f);
	};
	rep(i,0,4) rep(j,i+1,4) rep(k,j+1,4)
		mf(i, j, k, 6 - i - j - k);

	rep(i,4,sz(A)) {
		rep(j,0,sz(FS)) {
			F f = FS[j];
			if(f.q.dot(A[i]) > f.q.dot(A[f.a])) {
				E(a,b).rem(f.c);
				E(a,c).rem(f.b);
				E(b,c).rem(f.a);
        std::swap(FS[j--], FS.back());
				FS.pop_back();
			}
		}
		int nw = sz(FS);
		rep(j,0,nw) {
			F f = FS[j];
#define C(a, b, c) if (E(a,b).cnt() != 2) mf(f.a, f.b, i, f.c);
			C(a, b, c); C(a, c, b); C(b, c, a);
		}
	}
	for (F& it : FS) if ((A[it.b] - A[it.a]).cross(
		A[it.c] - A[it.a]).dot(it.q) <= 0) std::swap(it.c, it.b);
	return FS;
};

template<class P, class F>
void delaunay(vector<P>& ps, F trifun) {
	if (sz(ps) == 3) { int d = (ps[0].cross(ps[1], ps[2]) < 0);
		trifun(0,1+d,2-d); }
	vector<P3> p3;
	for (P p : ps) p3.emplace_back(p.x, p.y, p.dist2());
	if (sz(ps) > 3) for(auto t:hull3d(p3)) if ((p3[t.b]-p3[t.a]).
			cross(p3[t.c]-p3[t.a]).dot(P3(0,0,1)) < 0)
		trifun(t.a, t.c, t.b);
}

int main() {
  std::cin.tie(0)->sync_with_stdio(0);
  int n; std::cin >> n;
  std::vector<int> radii(n);
  vector<Point<int>> centers(n);
  for (int i = 0; i < n; ++i) {
    std::cin >> centers[i].x >> centers[i].y >> radii[i];
  }

  vector<std::tuple<double, int, int>> edges;
  auto add_edge = [&](int i, int j) {
    double d = (centers[i] - centers[j]).dist();
    edges.emplace_back(d, i, j);
  };

  if (n <= 1000) {
    for (int i = 0; i < n; ++i) {
      for (int j = i+1; j < n; ++j) {
        add_edge(i, j);
      }
    }
  } else {
    delaunay(centers, [&](int a, int b, int c) {
      add_edge(a, b);
      add_edge(b, c);
      add_edge(c, a);
    });
  }

  std::sort(edges.begin(), edges.end());

  auto uf = UFDS(n);
  std::cout << "YES\n";
  for (auto [d, i, j] : edges) {
    if (uf.unite(i, j)) {
      std::cout << centers[i].x << ' ' << centers[i].y << " " << centers[j].x << ' ' << centers[j].y << "\n";
    }
  }
}