代码求条(前方高能
查看原帖
代码求条(前方高能
322792
Alea楼主2023/10/25 09:20
#include <iostream>
#include <cmath>
#include <iomanip>
#include <set>
#include <algorithm>
#include <vector>
#define fs(X) fixed<<setprecision(X)
#define wf(X,Y) setw(X)<<setfill(Y)
using namespace std;
const double EPS=8e-5;
struct point{
    double x,y;
    point(double cx=0,double cy=0):x(cx),y(cy){}
};
point operator*(double a,point p){
    return point(p.x*a,p.y*a);
}
point operator+(point a,point b){
    return point(a.x+b.x,a.y+b.y);
}
point operator/(point a,double b){
    double x=a.x/b,y=a.y/b;
    if(fabs(a.x)<EPS||fabs(b)<EPS) x=0;
    if(fabs(a.y)<EPS||fabs(b)<EPS) y=0;
    return point(x,y);
}
struct nurbs{
    int n; // numbers of control points
    int k; // degree
    int m; // number of knots
    point P[25]; // control points
    double w[25]; // weight of points
    double t[25]; // knot vector
    double N(int i,int k,double u){ // Function N
        if(k==0) return (t[i]<=u&&u<t[i+1]?1.0:0.0);
        double co0,co1;
        if(fabs(t[i+k]-t[i])<EPS||fabs(u-t[i])<EPS) co0=0;
        else co0=(u-t[i])/(t[i+k]-t[i]);
        if(fabs(t[i+k+1]-u)<EPS||fabs(t[i+k+1]-t[i+1])<EPS) co1=0;
        else co1=(t[i+k+1]-u)/(t[i+k+1]-t[i+1]);
        return co0*N(i,k-1,u)+co1*N(i+1,k-1,u);
    }
    point C(double u){ // Function C (for a single curve)
        point num=point(0,0);
        double dem=0;
        for(int i=1;i<=n;i++){
            num=num+w[i]*N(i,k,u)*P[i];
            dem+=w[i]*N(i,k,u);
        }
        return num/dem;
    }
    void clear(){
        n=k=m=0;
        for(int i=0;i<25;i++) P[i].x=P[i].y=w[i]=t[i]=0;
    }
}curv[3];
set<pair<int,int>> pts[3];
void run(nurbs x,int o){
    for(double a=x.t[0];a<x.t[x.m];a+=EPS){
        point pt=x.C(a);
        pt.x*=1000,pt.y*=1000;
        pt.x=round(pt.x),pt.y=round(pt.y);
        pts[o].insert(make_pair(pt.x,pt.y));
    }
}
int main(){
    int T;
    cin>>T;
    for(int tid=1;tid<=T;tid++){
        pts[0].clear();
        for(int i=1;i<=2;i++) curv[i].clear(),pts[i].clear();
        for(int i=1;i<=2;i++){
            cin>>curv[i].n>>curv[i].m,curv[i].k=curv[i].m-curv[i].n-1;
            for(int j=1;j<=curv[i].n;j++) cin>>curv[i].P[j].x>>curv[i].P[j].y>>curv[i].w[j];
            for(int j=1;j<=curv[i].m;j++) cin>>curv[i].t[j];
        }
        for(int i=1;i<=2;i++) run(curv[i],i);
        set_intersection(pts[1].begin(),pts[1].end(),pts[2].begin(),pts[2].end(),inserter(pts[0],pts[0].begin()));
        vector<pair<int,int>> ans,shans;
        for(auto it:pts[0]) ans.push_back(it);
        shans.push_back(ans[0]);
        for(int i=1;i<ans.size();i++){
            if(abs(ans[i].first-ans[i-1].first)<5&&abs(ans[i].second-ans[i-1].second)<5) continue;
            shans.push_back(ans[i]);
        }
        cout<<"Case "<<tid<<": "<<shans.size()<<endl;
        for(auto it:shans) cout<<"("<<it.first/1000<<"."<<wf(3,'0')<<it.first%1000<<", "<<it.second/1000<<"."<<wf(3,'0')<<it.second%1000<<")"<<endl;
        cout<<endl;
    }
    return 0;
}
2023/10/25 09:20
加载中...