Haskell-Things / ImplicitCAD

A math-inspired CAD program in haskell. CSG, bevels, and shells; 2D & 3D geometry; 2D gcode generation...
https://implicitcad.org/
GNU Affero General Public License v3.0
1.39k stars 142 forks source link

Defining an implicit function #233

Open SmoothDragon opened 5 years ago

SmoothDragon commented 5 years ago

I'm not sure if this is possible, but it is certainly in the spirit of ImplicitCAD, so I would like to know if anyone knows how to do it. A Haskell answer is just as good as an extopenscad answer.

Suppose I want to define a new implicit function, in this case, the astroid. It is similar to the circle, but has a different L-norm.

astroid(x,y) = (x**(2/3)+y**(2/3))**(3/2) -1

This function is negative for (x,y) in the astroid region, zero on the perimeter and positive outside. I would now like to use this function in ImplicitCAD in this manner (replacing circle(1) with the new astroid shape):

This works:

linear_extrude(height(x,y) = 1 + cos(x)/2 + cos(y)/2 )
    circle(1);

I would like to make this:

linear_extrude(height(x,y) = 1 + cos(x)/2 + cos(y)/2 )
    astroid;
julialongtin commented 5 years ago

The code changes to do this are pretty small. I should make this a bit more convenient and generic, so you can do this from pure haskell, or pure extopenscad. I'm elbow deep in this code at the moment, so expect a more rigorous response in the next week or so. basically, you have to supply three functions: an argument parser, a bounding box calculator, and a function that given a coordinate, determines the implicit result of your object. I'll see if i can expose this conveniently.

for reference, see implicit/Primitives.hs, GetBox3.hs, and GetImplicit3.hs.

On Tue, Sep 24, 2019 at 4:12 PM SmoothDragon notifications@github.com wrote:

I'm not sure if this is possible, but it is certainly in the spirit of ImplicitCAD, so I would like to know if anyone knows how to do it. A Haskell answer is just as good as an extopenscad answer.

Suppose I want to define a new implicit function, in this case, the astroid. It is similar to the circle, but has a different L-norm.

astroid(x,y) = (x(2/3)+y(2/3))**(3/2) -1

This function is negative for (x,y) in the astroid region, zero on the perimeter and positive outside. I would now like to use this function in ImplicitCAD in this manner (replacing circle(1) with the new astroid shape):

This works:

linear_extrude(height(x,y) = 1 + cos(x)/2 + cos(y)/2 ) circle(1);

I would like to make this:

linear_extrude(height(x,y) = 1 + cos(x)/2 + cos(y)/2 ) astroid;

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/colah/ImplicitCAD/issues/233?email_source=notifications&email_token=AEAMAAQKKQCHXPM6CCLD6BLQLIU5ZA5CNFSM4I2BBJEKYY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4HNLDOZQ, or mute the thread https://github.com/notifications/unsubscribe-auth/AEAMAAXHWJ7J5EDXRCBL5STQLIU5ZANCNFSM4I2BBJEA .

SmoothDragon commented 5 years ago

I the implicit function approach to CAD is very interesting. If arbitrary equations can be used as the building blocks, a lot of things that are hard in OpenSCAD become easier. I'll look through the code you suggested.

SmoothDragon commented 4 years ago

I was able to finish my first hypocycloid project: https://github.com/SmoothDragon/HypocycloidVase. That would have been trivial with ImplicitCAD. I was dragging and rotating an xy curve along a spiral path. My next step to union such pieces seems to be creating all sorts of manifold problems for OpenSCAD. I'm looking forward to this feature when done. I'll take another crack at the source code in the mean time. I really should learn Haskell at some point.

canadaduane commented 4 years ago

This is really neat. I'd like to see more examples of this some day. I think it could lend itself well to stylistic layovers--imagine rough.js (it replaces basic shapes with "hand-drawn" styled shapes) but being able to swap in various styles in order to generate a 3D object.

julialongtin commented 4 years ago

I'm thinking i should add an implicit function passthrough to the scad engine.

EG:

raw3D(bbox=function_that_returns_an_r3, implicit(x,y,z)=function_that_returns_an_r); where x,y,z = position being evaluated.

p4l1ly commented 4 years ago

I don't know the openscad system but the problem is that real exponent of negative numbers is complex. You need to deal with it. Maybe you want to use abs like this? (abs x ** (2/3) + abs y ** (2/3))**(3/2) - 1

You told you're ok with a Haskell solution, so I've made a few of them:

import Graphics.Implicit

main = writeSTL 0.1 "foo.stl"
  $ extrudeRM
  0 -- no rounding
  (Left 0) -- constant twist
  (Left 1) -- constant scale
  (Left (0, 0)) -- constant translation
  (implicit astroid ((-1, -1), (1, 1))) -- the object, bounded by the box ((-1, -1), (1, 1))
  (Right$ \(x, y) -> 1 + cos(x)/2 + cos(y)/2)  -- non-constant height
  where
   astroid (x, y) = (abs x ** (2/3) + abs y ** (2/3))**(3/2) - 1

Or you can make a 3D astroid:

main = writeSTL 0.1 "foo.stl" $ implicit astroid3d ((-1, -1, -1), (1, 1, 1))
  where astroid3d (x, y, z) = (abs x**(2/3)+abs y**(2/3)+abs z**(2/3))**(3/2) -1

Or, if you want more flexibility, you may implement (and then modify for your needs) the extrusion very simply using implicit functions:

astro = writeSTL 0.1 "foo.stl" $ implicit astroidColumn ((-1, -1, 0), (1, 1, 10))
  where
    astroidColumn (x, y, z) = maximum
      [(abs x**(2/3)+abs y**(2/3))**(3/2) -1
      , -z
      , z - (1 + cos(x)/2 + cos(y)/2)
      ]

Or just write the 2D astroid

main = writeSVG 0.1 "foo.svg" $ scale (10, 10)$ implicit astroid ((-1, -1), (1, 1))
  where astroid (x, y) = (abs x**(2/3)+abs y**(2/3))**(3/2) -1

GLHF and try using Haskell, it gives you freedom. Btw, I've used my local fork of ImplicitCAD, which implements the change proposed here https://github.com/colah/ImplicitCAD/issues/234 . Not sure if the bug affects these examples.

SmoothDragon commented 4 years ago

Thanks for looking at this. Some exponents with rational coefficients, don't have complex answers. In this case, the power in question, (2/3), always returns a positive result. You can square your real number and then safely take the cube root of it. (Cube roots are well defined for any real number, so you could also take the cube first and then square it.) The outside power of (3/2) is well defined since you are summing positive real numbers. Thus the abs is unnecessary. It may be that Haskell calculates a real number first, and then exponentiates, in which case I might have to write it as two successive operations.

I will try playing around your Haskell code and seeing if I can get it to work. I like your idea.

p4l1ly commented 4 years ago

I know :) I just thought that abs x**(2/3) is more beautiful than (x^2)**(1/3). But it's a matter of taste.

p4l1ly commented 4 years ago

And it's not that Haskell calculates the cubic root first and then it calculates the power. If you use **, it calculates x**y as exp(y * ln(x)) so if x is negative, ln(x) fails, no matter what was y. It is a common behaviour in most programming languages (try out e.g. (-2)**(2/3) in Python3). On the other hand, if you want the power to be calculated as a product, you need to use ^: i.e. x^3 will be calculated as x*x*x (don't try this in Python, ^ means something else there).

balsoft commented 4 years ago

I'm not sure if I should open a separate issue, but after playing with implicit I've found some strange behaviour of bounding boxes.

Consider the following:

import Graphics.Implicit

hyperboloid (x, y, z) = x^2 + y^2 - z^2/5 -5

sheet (x, y, z) = z + 2

obj = implicit (hyperboloid * sheet) ((-10, -10, -2), (10, 10, 10))

main = writeSTL 0.5 "hyperboloid.stl" obj

I would expect to see a hyperboloid "standing" on a sheet.

However, I see the following:

hyperboloid_buggy

What I would expect would look something like

hyperboloid_expectation

And it only gets weirder.

Consider the following:

import Graphics.Implicit

hyperboloid (x, y, z) = x^2 + y^2 - z^2/5 -5

hyperboloidObj = implicit hyperboloid ((-10, -10, -2), (10, 10, 10))

sheetObj = rect3R 1 (-10, -10, -4) (10, 10, -4)

main = writeSTL 0.5 "hyperboloid.stl" (union [ hyperboloidObj, sheetObj ])

I would expect to see a hyperboloid "hovering" on top of a rectangle, or, at the very least, "standing" on it (taking the previous bug into account).

However, what I see is hyperboloid_buggy2

If I move the sheet into the bounding box of the hyperboloid, I see something very similar to the previous example.

Last example of strangeness:

import Graphics.Implicit

hyperboloid (x, y, z) = x^2 + y^2 - z^2/5 -5

sheet (x, y, z) = z + 5

hyperboloidObj = implicit hyperboloid ((-10, -10, -2), (10, 10, 10))

sheetObj = implicit sheet ((-10, -10, -5), (10, 10, 10))

main = writeSTL 0.5 "hyperboloid.stl" (union [ hyperboloidObj, sheetObj ])

I would expect to see a hyperboloid hovering on top of the sheet, instead I see

hyperboloid_buggy3

Note that there is a "hole" in the sheet, and that the hyperboloid always extends right to the sheet, however low I place it.

Changing the precision doesn't change much.

p4l1ly commented 4 years ago

@balsoft bounding box is only to advise the rendering engine, where to render: the engine can decide not to evaluate the function outside its bounding box to save some CPU time, but this is not guaranteed. Your function should not be negative outside the box, otherwise the behaviour is (let's say) undefined (depending on the current implementation of the engine).

Your sheet, nor hyperboloid are not positive outside their bounding boxes and that's the problem. Second, why are you multiplying the functions in the first example? If you want to make union, use min, if you want to make intersection, use max. Multiplication with plane only switches the signs below it, therefore you see the border anyways.

Here is something that works:

import Graphics.Implicit

hyperboloid (x, y, z) = x^2 + y^2 - z^2/5 -5
cropped_hyperboloid p@(x, y, z) = maximum [hyperboloid p, 2 - z, z - 10]

hyperboloidObj :: SymbolicObj3
hyperboloidObj = implicit cropped_hyperboloid ((-10, -10, -2), (10, 10, 10))

sheetObj = rect3R 1 (-10, -10, -5) (10, 10, -4)

main = writeSTL 0.5 "hyperboloid.stl"$ union [sheetObj, hyperboloidObj]
balsoft commented 4 years ago

@p4l1ly thank you very much! I think I have already figured out most of this, however, I'd prefer to have a clear error when the behavior is undefined, instead of rendering something that clearly doesn't make any sense.

I think the main issue here is that I didn't read the paper that implicit is based on, now I mostly understand the way it works and it's much easier to play around.