If someone is interested, here some m3d code:
function CalcDE(It3Dex: TPIteration3Dext; mctp: PMCTparameter): Double;
var dCT: TVec3D;
Rtmp, RD1, RD2, RD3, Rst: Double;
stmp, stmp2: Single;
ftmp, itt, mi: Integer;
CS, bIsC: LongBool;
begin
if mctp.IsCustomDE then //analytical DE method, based on one additional parameter
with mctp^ do
begin
if DEoption = 20 then It3Dex.RStopD := msDEstop * StepWidth * 1.03;
Result := mMandFunctionDE(@It3Dex.C1) * dDEscale;
end
else with It3Dex^ do
begin //3d gradient
CS := CalcSIT;
mctp.mMandFunction(@C1);
CalcSIT := False;
mi := MaxIt;
MaxIt := ItResultI;
RStop := RStop * 8;
Rtmp := Rout;
Rst := mctp.mctDEoffset;
mCopyVec(@dCT, @C1);
mAddVecWeight(@C1, @mctp.Vgrads[2], Rst);
mctp.mMandFunction(@C1);
RD1 := Sqr(Rtmp - Rout);
mCopyAddVecWeight(@C1, @dCT, @mctp.Vgrads[1], Rst);
mctp.mMandFunction(@C1);
RD2 := Sqr(Rtmp - Rout);
mCopyAddVecWeight(@C1, @dCT, @mctp.Vgrads[0], Rst);
mctp.mMandFunction(@C1);
RD3 := Sqr(Rtmp - Rout);
Result := Rtmp * Ln(Rtmp) * mctp.dDEscale / (Sqrt(RD1 + RD2 + RD3) + (Rst * 0.06));
Rout := Rtmp;
MaxIt := mi;
RStop := mctp.dRStop;
mCopyVec(@C1, @dCT);
CalcSIT := CS;
end;
if Result < mctp.msDEstop * 0.25 then Result := mctp.msDEstop * 0.25;
if mctp.FormulaType > 0 then //DEcombinate, second formula
begin
itt := It3Dex.ItResultI;
Rst := It3Dex.Otrap;
stmp2 := It3Dex.SmoothItD;
ftmp := mctp.FormulaType;
mctp.FormulaType := 0;
stmp := mctp.dDEscale;
mctp.dDEscale := mctp.dDEscale2;
bIsC := mctp.IsCustomDE;
mctp.IsCustomDE := mctp.IsCustomDE2;
It3Dex.nHybrid[0] := 999999;
Rtmp := calcDE(It3Dex, mctp);
case ftmp of
1: if Rtmp < Result then Result := Rtmp else begin It3Dex.ItResultI := itt; It3Dex.SmoothItD := stmp2; It3Dex.Otrap := Rst; end;
2: if Rtmp > Result then Result := Rtmp else begin It3Dex.ItResultI := itt; It3Dex.SmoothItD := stmp2; It3Dex.Otrap := Rst; end;
3: begin
Result := Exp((Ln(Result) + Ln(Rtmp)) * 0.5);
It3Dex.ItResultI := Min(itt, It3Dex.ItResultI);
It3Dex.SmoothItD := (It3Dex.SmoothItD + stmp2) * 0.5;
It3Dex.Otrap := (It3Dex.Otrap + Rst) * 0.5;
end;
end;
It3Dex.nHybrid[0] := 0;
mctp.FormulaType := ftmp;
mctp.dDEscale := stmp;
mctp.IsCustomDE := bIsC;
end;
end;
Raymarching part:
RSFmul := 1;
RLastStepWidth := dTmp * sZstepDiv; //dTmp is the distance estimation, here on the first calculation point
repeat
if Iteration3Dext.ItResultI >= iMaxIt2 then // if inside while stepping, step back
begin
dT1 := -0.5 * RLastStepWidth;
mZZ := mZZ + dT1;
mAddVecWeight(@Iteration3Dext.C1, @VgradsFOV, dT1);
msDEstop := DEstop * (1 + mZZ * mctDEstopFactor);
dTmp := CalcDE(@Iteration3Dext, @MCTparas);
RLastStepWidth := - dT1;
end;
if (Iteration3Dext.ItResultI < iMinIt) or
((Iteration3Dext.ItResultI < iMaxIt2) and (dTmp >= msDEstop)) then //##### next step ######
begin
RLastDE := dTmp;
dTmp := dTmp * sZstepDiv * RSFmul;
dT1 := MaxCS(msDEstop, 0.4) * mctMH04ZSD;
if dT1 < dTmp then
begin
if not bFirstStep then StepCount := StepCount + dT1 / dTmp else StepCount := StepCount + Random;
dTmp := dT1;
end
else StepCount := StepCount + 1;
if bFirstStep then
begin
seed := 214013 * seed + 2531011;
bFirstStep := False;
dTmp := ((seed shr 16) and $7FFF) * 0.000030517578125 * dTmp;
end;
RLastStepWidth := dTmp;
mZZ := mZZ + dTmp;
mAddVecWeight(@Iteration3Dext.C1, @VgradsFOV, dTmp);
msDEstop := DEstop * (1 + mZZ * mctDEstopFactor);
dTmp := CalcDE(@Iteration3Dext, @MCTparas);
if dTmp > RLastDE + RLastStepWidth then dTmp := RLastDE + RLastStepWidth; //helps a little on julias etc
if RLastDE > dTmp + s1em30 then //a little down regulation
begin
dT1 := RLastStepWidth / (RLastDE - dTmp);
if dT1 < 1 then
RSFmul := maxCS(0.5, dT1)
else
RSFmul := 1;
end
else RSFmul := 1;
end
...
until (mZZ > Zend) or PCalcThreadStats.pLBcalcStop^;