曲线之间的3D插值 [英] 3D interpolation between curves

查看:47
本文介绍了曲线之间的3D插值的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我有一组与温度有关的曲线.即曲线Mat1适用于温度310°C,Mat2适用于温度420°C.

如您所见,以对数刻度显示时,数据看起来更好;

现在,我需要通过内插Mat1和Mat2曲线来获得温度为370C的Mat3曲线.解决这个问题的最佳方法是什么?我猜想我可能需要进行某种3D插值.还需要考虑数据的性质(对数行为).

这是Mat1的数据

 <代码> 9.43E + 06 6.00E + 043.96E + 06 6.20E + 041.78E + 06 6.40E + 048.52E + 05 6.60E + 044.28E + 05 6.80E + 042.25E + 05 7.00E + 041.23E + 05 7.20E + 046.95E + 04 7.40E + 044.05E + 04 7.60E + 042.43E + 04 7.80E + 041.49E + 04 8.00E + 049.39E + 03 8.20E + 04 

这是Mat2的数据

  5.14E + 08 4.80E + 041.35E + 08 5.00E + 044.36E + 07 5.20E + 041.64E + 07 5.40E + 046.90E + 06 5.60E + 043.18E + 06 5.80E + 041.58E + 06 6.00E + 048.35E + 05 6.20E + 044.64E + 05 6.40E + 042.69E + 05 6.60E + 041.62E + 05 6.80E + 041.01E + 05 7.00E + 046.47E + 04 7.20E + 044.25E + 04 7.40E + 042.86E + 04 7.60E + 041.96E + 04 7.80E + 041.37E + 04 8.00E + 049735.23 8.20E + 04 

任何帮助将不胜感激.

我正在为另外两条曲线添加数据;

在21°C的温度下弯曲

 <代码> 3.98E + 07 6.30E + 041.58E + 07 6.40E + 044.03E + 06 6.60E + 041.47E + 06 6.80E + 046.57E + 05 7.00E + 043.37E + 05 7.20E + 041.91E + 05 7.40E + 041.16E + 05 7.60E + 047.49E + 04 7.80E + 045.04E + 04 8.00E + 043.52E + 04 8.20E + 042.53E + 04 8.40E + 041.87E + 04 8.60E + 041.41E + 04 8.80E + 041.08E + 04 9.00E + 048.47E + 03 9.20E + 04 

在537°C的温度下弯曲

 <代码> 7.91E + 06 3.80E + 043.29E + 06 4.00E + 041.51E + 06 4.20E + 047.48E + 05 4.40E + 043.95E + 05 4.60E + 042.20E + 05 4.80E + 041.28E + 05 5.00E + 047.77E + 04 5.20E + 044.87E + 04 5.40E + 043.14E + 04 5.60E + 042.08E + 04 5.80E + 041.41E + 04 6.00E + 049.73E + 03 6.20E + 046.85E + 03 6.40E + 04 

有关曲线的更多信息-这些是交变应力(y轴),材料在不同温度下的破坏循环次数(x轴).

谢谢.

解决方案

我设法获得了简单的示例.首先,必须对您的数据进行排序,以便必须按温度对测量进行排序,并且必须按 y (应力)对每个测量进行排序.我用升序.第一算法:

  1. 计算BBOX

    仅计算所有测量的最小和最大 x,y 坐标.这将用于对数刻度和线性刻度之间的转换以及对齐.

  2. 重新采样并对齐所有测量值

    因此将所有测量值转换为其样本具有相同的 y 值(跨所有测量值).我使用了统一采样的 y 轴.因此,简单的步骤是(ymax-ymin)/(n-1),其中 n 是重采样数据的点数.因此,所有度量都将具有相同的大小,并且所有 y 值在同一索引上的所有度量都将相同.丢失的 x 数据将填充为 0 .

    可以以线性比例进行重新采样.我使用了

    这里是此的C ++/VCL示例(只是忽略VCL内容):

     <代码>//$$ ----表格CPP ----//---------------------------------------------------------------------------#include< vcl.h>#include< math.h>#pragma hdrstop#include" win_main.h"//---------------------------------------------------------------------------#pragma软件包(smart_init)#pragma资源"* .dfm"TForm1 * Form1;//---------------------------------------------------------------------------int xs,ys;//屏幕分辨率图形:: TBitmap * bmp;//用于渲染的后缓冲区位图//---------------------------------------------------------------------------//从这里开始重要的事情//---------------------------------------------------------------------------float in [4] [40] =//输入测量格式为:{temperature,x0,y0,x1,y1 ...,-1}{{21.0,3.98E + 07,6.30E + 04,1.58E + 07,6.40E + 04,4.03E + 06,6.60E + 04,1.47E + 06,6.80E + 04,6.57E + 05,7.00E + 04,3.37E + 05,7.20E + 04,1.91E + 05,7.40E + 04,1.16E + 05,7.60E + 04,7.49E + 04、7.80E + 04,5.04E + 04,8.00E + 04,3.52E + 04、8.20E + 04,2.53E + 04、8.40E + 04,1.87E + 04,8.60E + 04,1.41E + 04、8.80E + 04,1.08E + 04、9.00E + 04,8.47E + 03,9.20E + 04,-1.0},{310.0,9.43E + 06,6.00E + 04,3.96E + 06,6.20E + 04,1.78E + 06,6.40E + 04,8.52E + 05,6.60E + 04,4.28E + 05,6.80E + 04,2.25E + 05,7.00E + 04,1.23E + 05,7.20E + 04,6.95E + 04,7.40E + 04,4.05E + 04、7.60E + 04,2.43E + 04、7.80E + 04,1.49E + 04,8.00E + 04,9.39E + 03,8.20E + 04,-1.0},{420.0,5.14E + 08,4.80E + 04,1.35E + 08,5.00E + 04,4.36E + 07,5.20E + 04,1.64E + 07,5.40E + 04,6.90E + 06,5.60E + 04,3.18E + 06,5.80E + 04,1.58E + 06,6.00E + 04,8.35E + 05,6.20E + 04,4.64E + 05,6.40E + 04,2.69E + 05,6.60E + 04,1.62E + 05,6.80E + 04,1.01E + 05,7.00E + 04,6.47E + 04,7.20E + 04,4.25E + 04、7.40E + 04,2.86E + 04,7.60E + 04,1.96E + 04、7.80E + 04,1.37E + 04,8.00E + 04,9735.23,8.20E + 04,-1.0},{537.0,7.91E + 06,3.80E + 04,3.29E + 06,4.00E + 04,1.51E + 06,4.20E + 04,7.48E + 05,4.40E + 04,3.95E + 05,4.60E + 04,2.20E + 05,4.80E + 04,1.28E + 05,5.00E + 04,7.77E + 04,5.20E + 04,4.87E + 04,5.40E + 04,3.14E + 04、5.60E + 04,2.08E + 04、5.80E + 04,1.41E + 04,6.00E + 04,9.73E + 03、6.20E + 04,6.85E + 03,6.40E + 04,-1.0}};//---------------------------------------------------------------------------//温度和输出数据//---------------------------------------------------------------------------const n = 40;//用指向重新映射曲线float dat [4] [2 + n + n];//重新采样输入曲线浮出[2 + n + n];//内插曲线浮动xmin,xmax,ymin,ymax;//BBOXvoid resample(float * out,float * in,float y0,float y1)//对y进行重新采样并将其对齐到范围和n个点并将其存储到out{浮点数t,d1,d2,a0,a1,a2,a3,x,y,x0,x1,x2,x3;int i,ii,i0,i1,i2,i3,nn;//扫描in []中有多少个点对于(nn = 0,i = 1; in [i]> = 0.0; i + = 2)nn ++;//将输入曲线重新采样到n个点out [0] = in [0];//复制Tout [n + n + 1] =-1;//数据结束对于(i = 0; i  = y)中断(i1 = 0; i1  y中)i1--;//近邻索引i0 = i1-1;如果(i0< 0)i0 = 0;i2 = i1 + 1;如果(i2> = nn)i2 = nn-1;i3 = i1 + 2;如果(i3> = nn)i3 = nn-1;//转换为数组索引i0 = 1 + i0 + i0;i1 = 1 + i1 + i1;i2 = 1 + i2 + i2;i3 = 1 + i3 + i3;//参数基于y值d1 = y-in [i1 + 1];d2 = in [i2 + 1] -in [i1 + 1];if(fabs(d2)> 1e-6)t = d1/d2;否则t = 0.0;//指向插值x0 = in [i0];x1 = in [i1];x2 = in [i2];x3 = in [i3];//x的三次插值d1 = 0.5 *(x2-x0);d2 = 0.5 *(x3-x1);a0 = x1;a1 = d1;a2 =(3.0 *(x2-x1))-(2.0 * d1)-d2;a3 = d1 + d2 +(2.0 *(-x2 + x1));x = a0 +(a1 * t)+(a2 * t * t)+(a3 * t * t * t);如果(x <0.0)x = 0.0;//只是为了确保不会弄乱数据//复制点out [ii + 0] = x;out [ii + 1] = y;}}//---------------------------------------------------------------------------void interpolate(float * out,float T)//根据与温度T匹配的dat [4] []作为n点曲线插入out []{//dat [] []必须按T,x,y升序int i,ii;//有效的T范围是< dat [1] [0],dat [2] [0]>浮点数t,d1,d2,a0,a1,a2,a3,x,x0,x1,x2,x3,t0,t1,t2,t3;out [0] = T;//复制Tout [n + n + 1] =-1;//数据结束//来自T的参数t =(T-dat [1] [0])/(dat [2] [0] -dat [1] [0]);t0 = dat [0] [0];t1 = dat [1] [0];t2 = dat [2] [0];t3 = dat [3] [0];//曲线之间的三次插值对于(i = 0; i  = 0.0;){x = dat [i];i ++;y = dat [i];i ++;如果(x <= 0.0)继续;如果(_reset){xmin = xmax = x;ymin = ymax = y;_reset = false;}如果(xmin> x)xmin = x;如果(xmax <x)xmax = x;如果(ymin> y)ymin = y;如果(ymax <y)ymax = y;}}//---------------------------------------------------------------------------void toscr(float& x,float& y)//将x,y从绘图数据转换为屏幕坐标(仅用于渲染){浮点数x0,dx,y1,dy;//范围< 0,1>//x =(x-xmin)/(xmax-xmin);//线性//y =(y-ymin)/(ymax-ymin);//线性(x> = xmin)≤x= log(x/xmin)/log(xmax/xmin):x = 0.0;//对数(y> = ymin)≤y= log(y/ymin)/log(ymax/ymin):y = 0.0;//对数//看法x0 = 0.1 * xs;dx = 0.8 * xs;y1 = 0.9 * ys;dy = 0.8 * ys;//[像素]x = x0 + x * dx;y = y1-y * dy;}//---------------------------------------------------------------------------void plot(float * dat,TColor col)//渲染测量数据(仅用于渲染){我浮点数x,y,r = 2;//曲线bmp->画布->笔-> Color = col;bmp->画布->字体-> Color = col;对于(e = 1,i = 1; dat [i]> = 0.0;){x = dat [i];i ++;y = dat [i];i ++;如果(x <= 0.0)继续;toscr(x,y);如果(e){bmp-> Canvas-> TextOutA(x,y,AnsiString().sprintf(%.0f C",dat [0]));bmp-> Canvas-> MoveTo(x,y);e = 0;}否则bmp-> Canvas-> LineTo(x,y);}//点对于(i = 1; dat [i]> = 0.0;){x = dat [i];i ++;y = dat [i];i ++;如果(x <= 0.0)继续;toscr(x,y);bmp-> Canvas-> Ellipse(x-r,y-r,x + r,y + r);}}//---------------------------------------------------------------------------void draw()//只渲染我的应用{bmp-> Canvas-> Brush-> Color = clWhite;bmp-> Canvas-> FillRect(TRect(0,0,xs,ys));plot(dat [0],clRed);plot(dat [1],clGreen);plot(dat [2],clBlue);plot(dat [3],clBlack);情节(out,clMaroon);表格1->画布->绘制(0,0,bmp);//bmp-> SaveToFile(" out.bmp");}//---------------------------------------------------------------------------__fastcall TForm1 :: TForm1(TComponent * Owner):TForm(Owner)//我的应用程序的初始化{//初始化后缓冲bmp = new Graphics :: TBitmap;bmp-> HandleType = bmDIB;bmp-> PixelFormat = pf32bit;//这里准备数据(重要)我对于(i = 0; i <4; i ++)minmax(in [i],i == 0);对于(i = 0; i <4; i ++)重采样(dat [i],in [i],ymin,ymax);//在这里为T = 370 [C]创建新数据插值(out,370.0);//并将其包含在BBOX中进行渲染minmax(out,false);}//---------------------------------------------------------------------------void __fastcall TForm1 :: FormDestroy(TObject * Sender)//不重要,只是我的应用程序的析构函数{删除bmp;}//---------------------------------------------------------------------------void __fastcall TForm1 :: FormResize(TObject * Sender)//不重要,只是调整大小事件{xs = ClientWidth;ys = ClientHeight;bmp-> Width = xs;bmp-> Height = ys;画();}//-------------------------------------------------------------------------void __fastcall TForm1 :: FormPaint(TObject * Sender)//不重要,只是重画事件{画();}//--------------------------------------------------------------------------- 

    有关如何使用它的信息,请参见函数 TForm1 :: TForm1(TComponent * Owner).

    但是物理有效性值得怀疑,您应该通过进行5次测量来测试这种插值是否会导致有效数据.使用4对5进行插值,并检查它们是否重叠.如果不重叠,则可能需要进行其他调整,例如增加插值多项式的度数,或者也将对数标度用于重采样等...

    I have a set of curves that are temperature dependent. I.e curve Mat1 is for temperature 310C and Mat2 is for temp 420C.

    As you can see, the data looks better when put in a logarithmic scale;

    Now I need to get Mat3 curve for temperature 370C by interpolating Mat1 and Mat2 curves. What is the best way to get around doing this? I'm guessing that I might need to do some sort of 3D interpolation. The nature of the data (logarithmic behavior) also needs to be considered.

    Here's the data for Mat1

    9.43E+06    6.00E+04
    3.96E+06    6.20E+04
    1.78E+06    6.40E+04
    8.52E+05    6.60E+04
    4.28E+05    6.80E+04
    2.25E+05    7.00E+04
    1.23E+05    7.20E+04
    6.95E+04    7.40E+04
    4.05E+04    7.60E+04
    2.43E+04    7.80E+04
    1.49E+04    8.00E+04
    9.39E+03    8.20E+04
    

    Here's the data for Mat2

    5.14E+08    4.80E+04
    1.35E+08    5.00E+04
    4.36E+07    5.20E+04
    1.64E+07    5.40E+04
    6.90E+06    5.60E+04
    3.18E+06    5.80E+04
    1.58E+06    6.00E+04
    8.35E+05    6.20E+04
    4.64E+05    6.40E+04
    2.69E+05    6.60E+04
    1.62E+05    6.80E+04
    1.01E+05    7.00E+04
    6.47E+04    7.20E+04
    4.25E+04    7.40E+04
    2.86E+04    7.60E+04
    1.96E+04    7.80E+04
    1.37E+04    8.00E+04
    9735.23     8.20E+04
    

    Any help would be appreciated.

    Edit: I'm adding data for two additional curves;

    Curve at temperature 21C

    3.98E+07    6.30E+04
    1.58E+07    6.40E+04
    4.03E+06    6.60E+04
    1.47E+06    6.80E+04
    6.57E+05    7.00E+04
    3.37E+05    7.20E+04
    1.91E+05    7.40E+04
    1.16E+05    7.60E+04
    7.49E+04    7.80E+04
    5.04E+04    8.00E+04
    3.52E+04    8.20E+04
    2.53E+04    8.40E+04
    1.87E+04    8.60E+04
    1.41E+04    8.80E+04
    1.08E+04    9.00E+04
    8.47E+03    9.20E+04
    

    Curve at temperature 537C

    7.91E+06    3.80E+04
    3.29E+06    4.00E+04
    1.51E+06    4.20E+04
    7.48E+05    4.40E+04
    3.95E+05    4.60E+04
    2.20E+05    4.80E+04
    1.28E+05    5.00E+04
    7.77E+04    5.20E+04
    4.87E+04    5.40E+04
    3.14E+04    5.60E+04
    2.08E+04    5.80E+04
    1.41E+04    6.00E+04
    9.73E+03    6.20E+04
    6.85E+03    6.40E+04
    

    More info about the curves - These are alternating stress (y axis), number of cycles to failure (x axis) curves for a material at different temperatures.

    Thanks.

    解决方案

    I managed to get simple example working. First of all your data must be ordered so the measurements must be sorted by temperature and each measurement must be ordered by y (stress). I used ascending order. First algorith:

    1. compute BBOX

      simply compute min and max x,y coordinates of all measurements together. This will be used for conversion between logarithmic and linear scale and also for aligning.

    2. resample and align all measurements

      so convert all of your measurements to form that it's samples are at the same y values (across all measurements). I used uniformly sampled y axis. So simply step is (ymax-ymin)/(n-1) where n is number of points of the resampled data. So all measurements will have the same size and all the y values will be the same across measurement on the same index. Missing x data will be filled with 0.

      The resampling can be done in linear scale. I used piecewise cubic interpolation.

    3. create new measurement for new temperature

      so simply create new measurement again containing n-points. The y value is the same as before (so just copy it from any of the aligned measurements) and then just take 1 point from each of the 4 measurements corresponding to the same point as we are processing and cubicaly interpolate its position. However this must be done in logarithmic scale!

      The valid range of temperature is between the 2nd and 3th measurement temperature.

    Here preview using your data and 370 C:

    And here C++/VCL example for this (just ignore the VCL stuff):

    //$$---- Form CPP ----
    //---------------------------------------------------------------------------
    #include <vcl.h>
    #include <math.h>
    #pragma hdrstop
    #include "win_main.h"
    //---------------------------------------------------------------------------
    #pragma package(smart_init)
    #pragma resource "*.dfm"
    TForm1 *Form1;
    //---------------------------------------------------------------------------
    int xs,ys;              // screen resolution
    Graphics::TBitmap *bmp; // back buffer bitmap for rendering
    //---------------------------------------------------------------------------
    // here starts the important stuff
    //---------------------------------------------------------------------------
    float in[4][40]=        // input measureements format is: { temperature,x0,y0,x1,y1...,-1 }
        {{ 21.0,
        3.98E+07,6.30E+04,
        1.58E+07,6.40E+04,
        4.03E+06,6.60E+04,
        1.47E+06,6.80E+04,
        6.57E+05,7.00E+04,
        3.37E+05,7.20E+04,
        1.91E+05,7.40E+04,
        1.16E+05,7.60E+04,
        7.49E+04,7.80E+04,
        5.04E+04,8.00E+04,
        3.52E+04,8.20E+04,
        2.53E+04,8.40E+04,
        1.87E+04,8.60E+04,
        1.41E+04,8.80E+04,
        1.08E+04,9.00E+04,
        8.47E+03,9.20E+04,
        -1.0 },
        { 310.0,
        9.43E+06,6.00E+04,
        3.96E+06,6.20E+04,
        1.78E+06,6.40E+04,
        8.52E+05,6.60E+04,
        4.28E+05,6.80E+04,
        2.25E+05,7.00E+04,
        1.23E+05,7.20E+04,
        6.95E+04,7.40E+04,
        4.05E+04,7.60E+04,
        2.43E+04,7.80E+04,
        1.49E+04,8.00E+04,
        9.39E+03,8.20E+04,
        -1.0 },
        { 420.0,
        5.14E+08,4.80E+04,
        1.35E+08,5.00E+04,
        4.36E+07,5.20E+04,
        1.64E+07,5.40E+04,
        6.90E+06,5.60E+04,
        3.18E+06,5.80E+04,
        1.58E+06,6.00E+04,
        8.35E+05,6.20E+04,
        4.64E+05,6.40E+04,
        2.69E+05,6.60E+04,
        1.62E+05,6.80E+04,
        1.01E+05,7.00E+04,
        6.47E+04,7.20E+04,
        4.25E+04,7.40E+04,
        2.86E+04,7.60E+04,
        1.96E+04,7.80E+04,
        1.37E+04,8.00E+04,
        9735.23 ,8.20E+04,
        -1.0 },
        { 537.0,
        7.91E+06,3.80E+04,
        3.29E+06,4.00E+04,
        1.51E+06,4.20E+04,
        7.48E+05,4.40E+04,
        3.95E+05,4.60E+04,
        2.20E+05,4.80E+04,
        1.28E+05,5.00E+04,
        7.77E+04,5.20E+04,
        4.87E+04,5.40E+04,
        3.14E+04,5.60E+04,
        2.08E+04,5.80E+04,
        1.41E+04,6.00E+04,
        9.73E+03,6.20E+04,
        6.85E+03,6.40E+04,
        -1.0 }};             
    //---------------------------------------------------------------------------
    // temp and output data
    //---------------------------------------------------------------------------
    const n=40;                         // points to resmaple curves with
    float dat[4][2+n+n];                // resampled input curves
    float out[2+n+n];                   // interpolated curve
    float xmin,xmax,ymin,ymax;          // BBOX
    void resample(float *out,float *in,float y0,float y1)   // resample and align y to range and n points and store it to out
        {
        float t,d1,d2,a0,a1,a2,a3,x,y,x0,x1,x2,x3;
        int i,ii,i0,i1,i2,i3,nn;
        // scan how many points in[] has
        for (nn=0,i=1;in[i]>=0.0;i+=2) nn++;
        // resample input curves to n points
        out[0]=in[0];           // copy T
        out[n+n+1]=-1;          // end of data
        for (i=0;i<n;i++)
            {
            // y uniformly distributed and aligned in the dat array
            y=y0+((y1-y0)*float(i)/float(n-1));
            ii=1+i +i ;
            // check if range present
            if ((y<in[1+1])||(y>in[1+nn-1+nn-1+1]))
                {
                out[ii+0]=0.0;
                out[ii+1]=y;
                continue;
                }
            // find i1 so in[i1] <= y < in[i1+1]
            // linear search, can be replaced with binary search
            for (i1=0;i1<nn;i1++) if (in[1+i1+i1+1]>=y) break;
            if (in[1+i1+i1+1]>y) i1--;
            // neigboring indexes
            i0=i1-1; if (i0<  0) i0=   0;
            i2=i1+1; if (i2>=nn) i2=nn-1;
            i3=i1+2; if (i3>=nn) i3=nn-1;
            // convert to array index
            i0=1+i0+i0;
            i1=1+i1+i1;
            i2=1+i2+i2;
            i3=1+i3+i3;
            // parameter is based on y value
            d1=y-in[i1+1];
            d2=in[i2+1]-in[i1+1];
            if (fabs(d2)>1e-6) t=d1/d2; else t=0.0;
            // points to interpolate
            x0=in[i0];
            x1=in[i1];
            x2=in[i2];
            x3=in[i3];
            // cubic interpoaltion of x
            d1=0.5*(x2-x0);
            d2=0.5*(x3-x1);
            a0=x1;
            a1=d1;
            a2=(3.0*(x2-x1))-(2.0*d1)-d2;
            a3=d1+d2+(2.0*(-x2+x1));
            x=a0+(a1*t)+(a2*t*t)+(a3*t*t*t);
            if (x<0.0) x=0.0;   // just to be sure data is not messed up
            // copy point
            out[ii+0]=x;
            out[ii+1]=y;
            }
        }
    //---------------------------------------------------------------------------
    void interpolate(float *out,float T) // interpolate out[] as n point curve from dat[4][] matching temperature T
        {                                // dat[][] must be ordered ascending by T,x,y
        int i,ii;                        // valid T range is <dat[1][0],dat[2][0]>
        float t,d1,d2,a0,a1,a2,a3,x,x0,x1,x2,x3,t0,t1,t2,t3;
        out[0]=T;               // copy T
        out[n+n+1]=-1;          // end of data
        // parameter from T
        t=(T-dat[1][0])/(dat[2][0]-dat[1][0]);
        t0=dat[0][0];
        t1=dat[1][0];
        t2=dat[2][0];
        t3=dat[3][0];
        // cubic interpolation between curves
        for (i=0;i<n;i++)
            {
            // points to interpolate
            ii=1+i+i;
            x0=dat[0][ii];
            x1=dat[1][ii];
            x2=dat[2][ii];
            x3=dat[3][ii];
            // logarithm scale
            (x0>=xmin)?x0=log(x0/xmin)/log(xmax/xmin):x0=0.0;
            (x1>=xmin)?x1=log(x1/xmin)/log(xmax/xmin):x1=0.0;
            (x2>=xmin)?x2=log(x2/xmin)/log(xmax/xmin):x2=0.0;
            (x3>=xmin)?x3=log(x3/xmin)/log(xmax/xmin):x3=0.0;
            out[ii+1]=dat[0][ii+1]; // copy y
            // too much missing data
            if ((x1<=0.0)||(x2<=0.0)){ out[ii+0]=0; continue; }
            // mirror missing data
            if (x0<=0.0) x0=x1-((x2-x1)*(t1-t0)/(t2-t1));
            if (x3<=0.0) x3=x2+((x2-x1)*(t3-t2)/(t2-t1));
            // interpolate x
            d1=0.5*(x2-x0);
            d2=0.5*(x3-x1);
            a0=x1;
            a1=d1;
            a2=(3.0*(x2-x1))-(2.0*d1)-d2;
            a3=d1+d2+(2.0*(-x2+x1));
            x=a0+(a1*t)+(a2*t*t)+(a3*t*t*t);
            if (x<0.0) x=0.0;   // just to be sure data is not messed up
            else x=exp(x*log(xmax/xmin))*xmin; // back to linear scale
            out[ii+0]=x;
            }
        }
    //---------------------------------------------------------------------------
    void minmax(float *dat,bool _reset) // compute BBOX of the curves
        {
        int i;
        float x,y;
        for (i=1;dat[i]>=0.0;)
            {
            x=dat[i]; i++;
            y=dat[i]; i++;
            if (x<=0.0) continue;
            if (_reset){ xmin=xmax=x; ymin=ymax=y; _reset=false; }
            if (xmin>x) xmin=x;
            if (xmax<x) xmax=x;
            if (ymin>y) ymin=y;
            if (ymax<y) ymax=y;
            }
        }
    //---------------------------------------------------------------------------
    void toscr(float &x,float &y)   // convert x,y from plot data to screen coordinates (just for rendering)
        {
        float x0,dx,y1,dy;
        // range <0,1>
    //  x=(x-xmin)/(xmax-xmin);                         // linear
    //  y=(y-ymin)/(ymax-ymin);                         // linear
        (x>=xmin)?x=log(x/xmin)/log(xmax/xmin):x=0.0;   // logarithmic
        (y>=ymin)?y=log(y/ymin)/log(ymax/ymin):y=0.0;   // logarithmic
        // view
        x0=0.1*xs; dx=0.8*xs;
        y1=0.9*ys; dy=0.8*ys;
        // [pixels]
        x=x0+x*dx;
        y=y1-y*dy;
        }
    //---------------------------------------------------------------------------
    void plot(float *dat,TColor col)// renders measurement data (just for rendering)
        {
        int i,e;
        float x,y,r=2;
        // curve
        bmp->Canvas->Pen->Color=col;
        bmp->Canvas->Font->Color=col;
        for (e=1,i=1;dat[i]>=0.0;)
            {
            x=dat[i]; i++;
            y=dat[i]; i++;
            if (x<=0.0) continue;
            toscr(x,y);
            if (e)
                {
                bmp->Canvas->TextOutA(x,y,AnsiString().sprintf("%.0f C",dat[0]));
                bmp->Canvas->MoveTo(x,y);
                e=0;
                }
            else bmp->Canvas->LineTo(x,y);
            }
        // points
        for (i=1;dat[i]>=0.0;)
            {
            x=dat[i]; i++;
            y=dat[i]; i++;
            if (x<=0.0) continue;
            toscr(x,y);
            bmp->Canvas->Ellipse(x-r,y-r,x+r,y+r);
            }
        }
    //---------------------------------------------------------------------------
    void draw() // just render of my App
        {
        bmp->Canvas->Brush->Color=clWhite;
        bmp->Canvas->FillRect(TRect(0,0,xs,ys));
    
        plot(dat[0],clRed);
        plot(dat[1],clGreen);
        plot(dat[2],clBlue);
        plot(dat[3],clBlack);
    
        plot(out,clMaroon);
    
        Form1->Canvas->Draw(0,0,bmp);
    //  bmp->SaveToFile("out.bmp");
        }
    //---------------------------------------------------------------------------
    __fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner) // init of my app
        {
        // init backbuffer
        bmp=new Graphics::TBitmap;
        bmp->HandleType=bmDIB;
        bmp->PixelFormat=pf32bit;
        // here prepare data (important)
        int i;
        for (i=0;i<4;i++) minmax(in[i],i==0);
        for (i=0;i<4;i++) resample(dat[i],in[i],ymin,ymax);
        // here create new data for T=370[C]
        interpolate(out,370.0);
        // and also include it to the BBOX for rendering
        minmax(out,false);
        }
    //---------------------------------------------------------------------------
    void __fastcall TForm1::FormDestroy(TObject *Sender) // not important just destructor of my App
        {
        delete bmp;
        }
    //---------------------------------------------------------------------------
    void __fastcall TForm1::FormResize(TObject *Sender) // not important just resize event
        {
        xs=ClientWidth;
        ys=ClientHeight;
        bmp->Width=xs;
        bmp->Height=ys;
        draw();
        }
    //-------------------------------------------------------------------------
    void __fastcall TForm1::FormPaint(TObject *Sender) // not important just repaint event
        {
        draw();
        }
    //---------------------------------------------------------------------------
    

    See function TForm1::TForm1(TComponent* Owner) on how to use this.

    However physical validity is questionable You should test if this kind of interpolation leads to valid data by having 5 measurements. Use 4 to interpolate the 5th and check if they overlap If not then this might need additional tweaking like increasing the interpolation polynomial degree, or use log scale also for resampling etc ...

    这篇关于曲线之间的3D插值的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆